Orbital and interlayer Skyrmions crystals in bilayer graphene 
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A graphene bilayer in a transverse magnetic field has a set of Landau levels with energies E = 
±.y/N (N + l)hw* where uj* is the effective cyclotron frequency and N = 0,1, 2, ... All Landau levels 
but N = are four times degenerate counting spin and valley degrees of freedom. The Landau level 
N = has an extra degeneracy due to the fact that orbitals n = and n = 1 both have zero kinetic 
energies. At integer filling factors, Coulomb interactions produce a set of broken-symmetry states 
with partial or full alignement in space of the valley and orbital pseudospins. These quantum Hall 
pseudo-ferromagnetic states support topological charged excitations in the form of orbital and valley 
Skyrmions. Away from integer fillings, these topological excitations can condense to form a rich 
variety of Skyrme crystals with interesting properties. We study in this paper different crystal phases 
that occur when an electric field is applied between the layers. We show that orbital Skyrmions, 
in analogy with spin Skyrmions, have a texture of electrical dipoles that can be controlled by an 
in-plane electric field. Moreover, the modulation of electronic density in the crystalline phases are 
experimentally accessible through a measurement of their local density of states 

PACS numbers: 73.21.-b,73.22.Gk,78.70.Gq 



I. INTRODUCTION 

Bilayer graphene is a system consisting of two layers 
of graphene separated by a distance d = 3.337 A. In the 
Bernal stacking structure, one of the two honeycomb sub- 
lattice sites in each layer has a near neighbor in the other 
layer and one does not. In a transverse magnetic field, 
the two-dimensional electron gas (2DEG) develops a set 
of Landau levels with energies E° = ±y/ N (N + Tjtfcj* 
where N = 0,1,2,... and uj* = eB/m*c is the effec- 
tive cyclotron frequency. The effective mass given by 
to* = 2fi 2 7 1 /37QdQ where a is the lattice constant of 
graphene and 70 and 71 are in-plane nearest-neighbor 
and inter-plane hopping parameters. By comparison, the 
effective mass is zero in graphene and the Landau levels 
energies are then given by E° — ±^/2HvfVN/£ where 
I 2 = hc/eB is the magnetic length and vf — V3o,ojo/^fi- 
is the Fermi velocity. 

In the absence of an electric held between the lay- 
ers and Zeeman coupling, the Landau level N = in a 
graphene bilayer contains 8 states for each guiding center 
orbital. The extra degeneracy is due to the fact that or- 
bitals (we define this term in the next section) n — and 
n = 1 both have the same kinetic energy E — 0. Con- 
sequently, an electron in N = must be described by 
its spin, valley (or layer), and orbital quantum numbers 
in addition to its guiding center index X in the Landau 
gauge. When Coulomb interaction is considered, this ex- 
tra degeneracy produces a rich phase diagram for the 
bilayer graphene's 2DEG. More so, in fact, than in a 
semiconductor 2DEG. In a series of related papers^— , 
we have shown that the octet degeneracy is lifted by 
the Coulomb interaction. The broken-symmetry ground 
states that emerge can be described as quantum Hall 



pseudo-ferromagnets in the pseudospin language where 
fictitious spins are associated with the valley and or- 
bital indices. These new states have interesting trans- 
port properties such as an intra-Landau-level cyclotron 
mode and a layer pseudospin with a quadratic (uj ~ q 2 ) 
dispersion implying a vanishing superfluid density. 

It is well known that a quantum Hall ferromagnet 
(QHF) in a usual semiconductor 2DEG has topological 
excitations named spin Skyrmions*. A single Skyrmion 
spin texture has its spins aligned with the Zeeman held at 
infinity, reversed at the center of the Skyrmion, and has 
non zero XY spin components at intermediate distance 
which have a vortex-like configuration. Skyrmions carry 
electric charge. Calculations have shown that Skyrmion- 
anti-Skyrmion pairs have lower energy than electron-hole 
quasiparticles near filling factor v = 1 and dominate the 
transport properties of the QHF—. A quantum Hall bi- 
layer has, in addition to spin Skyrmions, topological exci- 
tations named pseudospin Skyrmions where the fictitious 
spin is associated with the layer index. We show in this 
paper that, in a graphene bilayer, there is a third pos- 
sibility: that of an orbital-pseudospin Skyrmion. This 
quasiparticle has an associated electric dipole texture in 
the plane of the layers and can be seen as the analog of 
a spin Skyrmion who carries a magnetic texture. 

We present several crystal states with pseudospin tex- 
ture that occur near integer filling factors in Landau level 
N = 0. We assume full spin polarization of the 2DEG and 
concentrate on Skyrme crystals with valley and/or orbital 
pseudospin textures. We allow for the presence of an elec- 
tric field between the layers that creates a charge imbal- 
ance and we study the evolution of the Skyrme crystals 
as a function of this electrical "bias" for different filling 
factors. Our goal is not to study the full phase diagram 
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of the bilayer but to focus on a small number of inter- 
esting crystal phases that are likely to occur near inte- 
ger filling factors. Our calculation shows that Skyrmions 
tend to crystallize in pairs and that an orbital-pseudospin 
texture is favored over valley-pseudospin texture at fill- 
ing factors v = —3,— 1,1,3. Valley-pseudospin Skyrmion 
crystals occur at filling factor v — —2,2 and involve a 
texture in both the n = and n = 1 orbitals. This possi- 
bility was discussed before, for an isolated Skyrmion, in 
Ref. [6J. We also show that Skyrmions with different elec- 
tric charge q can be distinguished on the basis of their 
density of states. Moreover, the real-space density pat- 
tern of a Skyrmion crystal is accessible by a measurement 
of its local density of states. 

This paper is organized as follows. In Sec. II, we 
present the effective two-band tight-binding model that 
we use to describe the graphene bilayer. In Sec. Ill, we 
derive the Hamiltonian of the bilayer graphene's 2DEG 
in the Hartree-Fock approximation. Sec. IV presents the 
pseudospin language used to describe the various crys- 
tal phases. In Sec. V, we introduce spin and orbital- 
pseudospin Skyrmions. We then present and discuss 
various Skyrmion crystal phases at filling factors v — 
—3,-2,-1 (or equivalently v — 1,2,3) in Sec. VI. The 
total and local densities of states are defined in Sec. VII 
and calculated for some of the crystal phases. The elec- 
tric dipole texture associated with an orbital Skyrmion 
crystal is computed in Sec. VIII. We conclude in Sec. IX 
with a discussion of some of the terms neglected in our 
simple tight-binding model. 



II. EFFECTIVE HAMILTONIAN 

We consider the graphene bilayer in the Bernal stack- 
ing arrangement^ represented in Fig. [TJ We denote the 
two basis atoms of the top layer by A\ and B\ and those 
of the bottom layer by A-i and B^ with atoms A\ sit- 
uated directly above atoms Bi- The bilayer is placed 
in an external transverse electric field in order to con- 
trol the electrical potential difference (i.e. the "bias") 
Ab between the layers that causes the charge imbalance. 
To simplify our analysis, we assume complete spin polar- 
ization of the electron gas and neglect trigonal warping 
(the 73 hopping in Fig. [I}. We also use an effective two- 
band model^ to describe the low-energy excitations of 
the bilayer in a quantizing magnetic field in the valleys 
K = (-47r/3a ,0) and K' = (47r/3a ,0). Although we 
will not consider these terms in the bulk of this paper, 
we could generalize this model by including the 74 hop- 
ping term as well as an additional term A representing 
the difference in the crystal field experienced by the in- 
equivalent atoms A and B in the same plane. With these 
approximations, we get the Hamiltonian: 



+ (A) As + Ci)aa t 
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FIG. 1: (Color online) Crystal structure and definition of the 
hopping parameters for the graphene bilayer. 

in the basis (A2,Bi) . In Eq. (p}, a,a^ are the ladder 
operators for the Landau levels and we have defined the 
parameters 
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are unitless constants and sgn denotes the signum func- 
tion. The effective cyclotron frequency is defined by 
lu* = eB/m*c with the effective electronic mass given 
by to* = 2fi 2 7i/37gaQ = 0.054too where mo is the 
bare electronic mass and ao = 2.46 A is the lattice pa- 
rameter of graphene. Note that, in the basis (A2,B\), 

Hh = {H^- 

In the case where 74 = A = As=0, the Landau level 
energies are given by 



E° = ±y/N (N + l)hu*, 



(6) 



with N = 0,1, 2, 3, ... All Landau levels are four time de- 
generate (including spin and valley degrees of freedom) 
with the exception of N = that is eight times degener- 
ate. With finite 74, A, Ab, we find for the spin up states 
of N = the following spinors and energies: 
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for the K valley and 



h o,x (r) i p0 

I ' ■ c, iT',0,X 

(r) \ F n 



(9) 

/3 A B + Ci,(io) 



for the K' valley. Note that we have neglected the 
Zeeman coupling since we assume complete spin po- 
larization and thus discard the spin degree of freedom 
in the rest of our analysis. The functions h n> x (r) = 
e -iXy/e ^ x _ I -JT^y are the eigenstates in the Lan- 
dau gauge A = (0,Bx, 0) with guiding center X and 
ip n (x) is the wave function of the one-dimensional har- 
monic oscillator. The magnetic length is given by £ — 
yJhcjeB = 256 /V~B A. We see from Eqs. 08) that, in 
addition to the spin and valley quantum numbers, there 
is in N = an extra degeneracy due to the fact that wave 
functions hg^x (r) and hi t x (r) both have zero kinetic en- 
ergy if As = 0. Throughout this paper, we will refer to 
these states as orbitals n = 0, 1 and use the symbol N 
for the Landau level index. 

From Eqs. (I7|)- ([TU1) . we see that, in the Landau level 
N = 0, the electrons are localized on the atoms Ai (bot- 
tom layer) in the K 1 valley and on the atoms B\ (top 
layer) in the K valley. In N — 0, the layer index is equiv- 
alent to the valley index. An external electric field lifts 
both the valley and the orbital degeneracies. The orbital 
degeneracy is lifted by the small corrections (3qAb and 
Ci as shown in Fig. [2] The values of the intra and inter- 
layer hoppings are given by 70 = 3.12 eV and 71 = 0.39 
eV. The other hopping terms as well as A are not so well 
known. It is difficult to get the relative signs of these 
terms from the litterature. Recent measurements of these 
parameters for bilayer graphene give 74 = 0.04 — 0.07 and 
A = 0.005 — 0.008 in units of the in-plane hopping 70. 
These values are discussed and referenced in Ref. 0- Tak- 
ing the minimal values for 74 and A, we have fta = 8. 
86 x 10^ 3 B,/3 4 = 1.31 X 10~ 5 B (with the magnetic field 
in Tesla) and A = 0.0156 eV so that Ci = 4. 042 x 10~ 4 B 
eV. 

We showed in Refs. that when £1 = 0, the phase 
diagram of the 2DEG at integer filling factors v € [—3, 4] 
contains phases with interlayer and/or inter-orbital co- 
herences. Because of the small interlayer spacing (d = 
3.337 A) in a graphene bilayer, interlayer coherence is 
rapidly lost when Ab increases i.e. for Ab/ (e 2 / ni) > 
0.001 according to our numerical calculations (k is the ef- 
fective dielectric constant at the position of the graphene 
layers). Above this value, inter-orbital coherence sets in 
when E° K0 X > E° K1X . From Eqs. PITU]). this is only 
possible at v = —1,3. Indeed, our calculations show that 
the phase diagram for v = —3, 1 has no orbital-coherent 
phase (when C,\ = 0) if the band parameters we use are 
correct. 

The precise values of the bias for the transitions be- 
tween the different liquid phases at integer filling fac- 
tors are very sensitive to the exact values of the hop- 
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FIG. 2: Non-interacting energy levels with spin up in Landau 
level N = 0. Note that when (1 = 0, level K, 1 is below level 
K, in energy. 



ping parameters. The same is true for the boundaries 
between various crystal phases at non-integer filling fac- 
tors. Moreover, the number of possible crystal phases is 
much larger than the number of possible liquid phases 
when one considers the various crystal lattices and the 
possibility of having more than one electron per unit cell 
with interlayer and/or orbital pseudospin textures. For 
this reason, we focus, in this paper, on the analysis of a 
few crystal phases with orbital or interlayer texture which 
are likely to appear in the phase diagram of the 2DEG 
in some range of values of £1 . We assume Ci = for all 
our calculations and discuss in the conclusion how the 
phase diagram is likely to be changed when £1 7^ 0. In 
our opinion reliable determination of the phase boundary 
characterizing the many possible crystalline phases will 
require experimental input. 

Note that in the absence of Landau level mixing and 
when maximal spin polarization is assumed, the ground 
states at filling factors v = —3,-2,-1 are equivalent to 
those at filling factors v — 1,2,3. It is thus sufficient for 
us to study the first three states v — —3,-2,-1. If the 
approximation of maximal spin polarization is not made, 
phases with reduced polarization become possible as the 
bias is increased. (The energy of half the spin down (up) 
states decreases (increases) with bias and levels cross- 
ing do occur). The ground states at v — —3, —2,-1 arc 
no longer equivalent to those at v — 1,2,3. The phase 
diagram is much more complex in this case but our cal- 
culations show that phases with orbital or interlayer co- 
herences are still present. 
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III. HARTREE-FOCK DESCRIPTION OF THE 
CRYSTAL PHASES 

We now add the Coulomb interaction to the non- 
interacting Hamiltonian. We assume that the magnetic 
field is strong enough so that we can neglect Landau level 
mixing. The Hartree-Fock Hamiltonian for the 2DEG in 
Landau level N = is then given by 



H 



HF 



(11) 



+^E 

a,b n 



/ j ni,n 2 ,r 



(q) 



*(p a ntn 2 (-q)Ktu (q) 

~N(p y] E^ -^ni,n4,n 3 ,n 2 

a,b m,...,m q 

x{p a n tn 2 (-q)>/4 a ,» 4 (q). 



(q) 



where iV^ = S/2ir£ 2 is the Landau level degeneracy (S 1 
is the 2DEG area) and all energies are now measured in 
units of e 2 / k£. The single-particle energies E a>n include 
capacitive contributions and are defined by 



E a ,n = T^Ab - af3 A B n 



v d „ d 
v„ — 

2£ £ 



+ Ci, (12) 



with a,b = ±1 the valley (or equivalently layer) index 
and n = 0, 1 the orbital index. (Our convention is that 
a = 1(— 1) for the K(K') valley). Because we work in 
Landau level N = only, we define ?=i/ + 46 [0, 8] 
as the number of filled levels in N — 0. In deriving 
Eq. (jlip . we have taken into account a neutralizing pos- 
itive background and put the capacitive energy in the 
third term on the right-hand side of Eq. (fT2|) . It follows 
that the q = contribution is absent in the Hartree term 
of Eq. (jlip . This convention is indicated by the bar over 
the summation. 

In Eq. (|TT|) . the density operator is defined by 



ntn 2 (^ = ^ E e-**<*+*> (13) 

v x u x 2 

XC l.X 1 .n 1 C b,X 2 ,n 2 Sx 1 ,X 2 +q y P, 



where c\ x creates an electron in the state (a, X, n) in 
the Landau gauge. The intralayer (H,X = H a > a 1 X a ' a ) 

and intcrlaycr (h,X = H a ^ b ,X a ^ Hartree and Fock 
interactions are defined by 



H nu n 2 ,n 3 ,n 4 (q) = -£#ni,n 3 (q) #n 3) n 4 (~q) ) (14) 



■^ni ,n 2 



(q) 



dpf 



Hni 



"2,13,14 



(P)e 



iq X pi ' 



, (15) 



and 
X 



14 (q) — ffni,Ti2,i3,i4 (q) e 

dp£ 2 - 



ni,n2,n 3 ,n 4 



(q) 



2^ - J ™1,«2,«3,™4 



(16) 
(p)e"^U7) 



where d — 3.337 A is the separation between the two 
graphene layers of the bilayer. The form factors which 
appear here, 



#0,0 (q) = exp 



(18) 
(19) 



#i,o (q) = 7= exp — - — , (20) 



#1,1 (q) = exp I ; — I I 1 



V2 

#o,i (q) = ^ *S= J exp I - 1 , (21) 

capture the character of the two different orbital states. 
Detailed expressions for the Hartree and Fock interac- 
tions can be found in Appendix A of Ref. [j§. 

The average values of the density operators 
(Pnfna (q)) are found by solving the Hartree-Fock 
equation of motion for the matrix Green's function 



G 



ntn 2 (q,r) = ^ E e- l M*+*') 



(22) 



x$x,x>-q y e 2 Glf n2 (X,X',t) 



with 



G# )Ba (X, X\ t) = - (Tc , x , m (r) c^, ^ (0)) . (23) 

When t = 0-, 

G^„ 2 ( q ,r = 0-) = (^ ni (q)). (24) 

The Hartre-Fock equation of motion for this single par- 
ticle Green's function is given by 



[Hiu> n - (E a<n - fi)] G a n n , (q,w„) = h8 Cifi S n . n '5 a . b (25) 
+ E EC "4, q - q') e - lqXq '^ /2 G- fc „, (q> n ) 

c,ri4 q' 

-EE C (n, "4, q - q') e-^^/'G^, (q', Wn ) , 

c,n4 q' 

with the Hartree and Fock potentials 

U" a ( n ! n 4, q) = E Hc >« ( n i' n 2, n, n 4 ; -q) (p^ c ,„ 2 (q)) , 

"1 ,"2 

(26) 

C ( n > n 4, q) = E Xc < a ( n i' n 4, n, n 2 ; -q) (p%" n3 (q)) • 

(27) 
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In these equations, H c=a = H, H c ^ a = H and similarly 
for X c>a . 

Equation ()25|) constitutes a set of self-consistent equa- 
tions that can be solved numerically using the procedure 
described in Ref. Oil We search amongst the many so- 
lutions of this equation for the one that minimizes the 
Hartree-Fock energy per electron 



E 



HF 



1 



No v 



= ^V^,„(p^(o)> 



(28) 



1 



a,b m,...,n4 q 

x(/C„ a (-«!)> (q)> 



2V ^ ^ ^ 

a,b nx,...,ri4 q 



^ Q: 6 (ni,n 4 ,n 3 ,n 2 ;q) 



x «,„ 2 (-q)) (ftt„ 4 (q)> , 
where Nq is the number of electrons in N = 0. 



IV. ORDER PARAMETERS AND PSEUDOSPIN 
DESCRIPTION 

The set of parameters {(Pn'in 2 (l)/} mn Y character- 
izes a particular ground state. In the uniform states 
studied in Refs. lU-El these parameters were nonzero 
for q = only. In this paper, however, we study crystal 
states occurring at non-integer filling factor v where we 
expect a finite fraction of the electrons, usually v — \v\ , 
to crystallize. The set of parameters {(Pn'i,n 2 (<3))} are 
then nonzero for q = G where G is a reciprocal lattice 
vector of the crystal lattice considered. 

Generally speaking, crystalline states should occur uni- 
versally near integer filling factors in order to maxi- 
mize the correlations among the lowest-energy elemen- 
tary charged excitations. When the charged objects be- 
come more dense at larger departures from integer filling 
factor they will begin to overlap. Eventually the various 
exotic crystalline states will quantum melt and the elec- 
trons will form a fluid state. The methods employed in 
this paper are not able to predict the stability range of 
the crystalline states. 

In Fourier space, the real electronic density in valley 
(or layer) a is given by 



(n a (G)) 



N v K n , m (-G) { P a n % (G)) . (29) 



We will refer to the inverse Fourier transform n a (r) of 
(n a (G)) as the density in the real space representation 
(RSR). We also make use of another expression for the 
density: 



(n„ (G)) = £ (pS£ (G)) 



(30) 



We refer to the inverse Fourier transform n a (r) of 
(n a (G)) as the density in the guiding- center representa- 
tion (GCR). By definition, (n a (G = 0)) = v a is just the 
filling factor in valley a. The form factors K n ^ m (G) are 
not taken into account in the GCR so that the character 
of the different orbitals n — 0, 1 is lost in the correspond- 
ing density. 

Interlayer coherence implies that (pn' b 7 f a (G)) ^ 
while inter-orbital coherence implies that 
Pnm=£n (G)^ 7^ 0. In the most general case, both 
interlayer and inter-orbital coherences are present and 
^Pn r$n (G)^ 7^ 0. The different phases are best de- 
scribed by using a pseudospin language. For the orbital 
pseudospin, S, we associate the up state with the n = 
orbital and the down state with the n = 1 orbital so 
that in valley a : 



Sa, Z (G) 



1 



[(Po',0 (G)) — (pi'i (G))] , (31) 



S a ,_i_(G) = Sa^x + S^yf , (32) 

Sa,+ (G) - §a,x + iSa,y = (Po'l (G)} ■ (33) 

For the interlayer pseudospin, P, we associate the up 
state with the K layer and the down state with the K' 
layer so that for orbital n, we have 



1 



«f (G)>- (p 



Pn,z (G) 

Px,„(G) = P„, x x + P„,3,y , 

Pn,-\- (G) — Pn,x H~ i>Pn,y — 



(G) 



Pnn ( G ) 



(34) 
(35) 
(36) 



Note that these fields are defined in the GCR. To get 
them in the RSR, we multiply each (p%m (G)) in these 
definitions with N v K n _ m (— G). From now on, we use 
the notation S, P, n to refer to the fields in the GCR and 
the notation S,P,n for the fields in the RSR. The two 
views give separate interesting insights into the nature of 
the crystal states, but the RSR is more closely related to 
experimental probes like a scanning tunneling microscope 
(STM). 

An exact description of the state of one electron in 
N = is given by the four complex components of the 

spinor (c^ K X Q , c^ x l , J K , xfi , c^, x l J . Note that, in or- 
der to limit the range of possible states, we will restrict 
our attention to circumstances in which the N — states 
are maximally polarized. As we mentioned in Sec. II, 
we expect that partially polarized states will be com- 
mon at large interlayer potentials. The ideas explained 
here are readily generalized to include this possibility. 
For a CP 3 spinor, the norm and the absolute phase are 
fixed so that a given electronic state is defined by 6 
independent components or angles 11 . The 12 classical 
fields P„= ,i (G) , S±x (G) that we defined in Eqs. (JUJ- 
I36p have a simple physical interpretation but they do not 
provide a full description of a given phase. Moreover, 
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these fields are not independent variables and their norm 
is not fixed. Both the modulus and the orientation of 
these pseudospins may vary in space. In fact, the sum 
rule 



EEEI«™( G )>I 



(37) 



a,b n.m G 



that applies when the many-electron state is approxi- 
mated by a single Slater determinant becomes, in pseu- 
dospin language, 



E 



G 



\\ PK {G) +PK , (G)| 2 



(38) 



+2 
+2 
+2 



P z , (G) - P z<1 (-G) 

2 



S K (G) 
Po (G) 



K,K' frt s 



S K > (G) 
Pi(G 



K,K' fr ,< 

Pub ( G , 



, where we have defined p a = Pq'q- 

Skyrmion crystals with intervalley pseudospin textures 
have been studied extensively in semiconductor 2DEG as 
well as in graphene monolayers^—. In bilayer graphene, 
we have the additional possibility of orbital pseudospin 
texture. This type of texture is particularly interesting 
because it gives rise to textures of electric dipoles in the 
plane of the layer. As shown in Refs. 013 the cou- 
pling of the 2DEG with an electric field E (r) = — V(j> (r) 
in the plane of the layers can be written as H ext = 
— e J dr n (r) <j> (r) where n (r) is the Fourier transform of 
the total density n K (G) + n K > (G) (see Eq. ([55)1. With 
the form factor defined in Eqs. (|18I21[) and j = K,K' , 



eN, 



S 



EE 



3 G 

G 2 £ 2 



G 2 



Pi (-G) 



(39) 



Pi,, (- G ) 



- V2t (GJ p j>x (-G) - G y £ p j>v (-G))] <f> (G) , 

where we have defined p - (G) = exp (— G 2 £ 2 /A) pj (G). 
In real space, 



H ext = -eN v E / dr \Pi ( r ) ^ ( r ) 



(p j (r)£ 2 -2p jjZ (r)^)(V-E(r)) 



(40) 



+V2£eN v / dr [p j}X (r) E x (r) - p^ y (r) E y (r)] , 



so that we can identify 

d a (G) = -eV^V^ 2 / 4 (41) 
x«p 0)ie (G))x-<p 0)J ,(G))y) 

with the Fourier transform of an electric dipole field in 
layer a. The orientation of the dipole vector at each point 
in space is simply related to the orientation of the orbital 
pseudospin vector. It follows that crystals with orbital 
pseudospin textures then have electric-dipole textures. 
Orbital Skyrmion crystals are the electric analog of spin 
Skyrmions crystals in which it is the magnetization that 
varies in space. Note that in a uniform electric field, 
the second term in Eq. (140f) is zero and H ex t gives the 
monopole and dipole terms of the interaction energy of 
the electrons with the external electric field. 



V. ISOLATED SKYRMIONS 

Before we can analyze the results of the numerical 
calculations for the Skyrmion crystal states, we need 
to know the density and pseudospin patterns associated 
with a single orbital Skyrmion located at r = 0. We use, 
in this section, the symmetric gauge which is more conve- 
nient for this problem. We take A = (— By/2, Bx/2, 0) 
for the vector potential. The eigenfunctions of the ki- 
netic Hamiltonian H = (p + eA/c) 2 /2mo (where —e is 
the charge of an electron and mo the electronic mass) are 
given by 

1 



h n =o,m (r) — 



-irmp —r /it 



V2ir2 m m\£ 

with m = 0, 1,2, ... for Landau level n = and by 

1 1 



(42) 



h n =l,m (r) 



(43) 



xe 



-imtp p -r 2 1 'it 2 jA m \ 



i+0 



2£ 2 J 



with to = —1, 0, 1, 2, ...for Landau level n = 1 and L™ (x) 
is a generalized Laguerre polynomial. (Note that h n , m 
is a different function than h n x introduced previously). 
Figure [3] shows the density profile n\ (r) = \h\. m (r)| 2 for 
the eigenstates with to = —1,0,1. 

The index to > —n gives the angular momentum i.e. 

L z h n . m (r) = -hmh n ^ m (r) , (44) 

while the energy of each state (n, to) is given by 

E° n = (n + 1/2) fkj c , (45) 

where u> c is the cyclotron frequency. Note that for a filled 
level, we have 



(r) 
m (r) 



oc 

E iv 

m=0 



1 



E 



\ h i,m (r)|' 



2tt£ 2 ' 
1 



2ir£ 2 ' 



(46) 
(47) 
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FIG. 3: (Color online) Density profile for the wave functions 



A. Spin Skyrmion 

In a semiconductor 2DEG, the state \S) corresponding 
to the addition of one spin Skyrmion with topological 
charge Q = 1 at r = to the spin-polarized ground state 

oo 

\ gs ) = n 4n,t i°) at v = 1 ( Landau ievei N = ° is 



implicitly assumed) is given by 

oo 

p>=n[ 



UmC \,m+l + VmC \,m 



m=0 

with the constraint that 

\Um\ 2 



io lo> 



The ant i- Skyrmion state is given by 

oo 

1-45) = J] [u m cl m + v m c\ m+1 ] |0) 



(48) 



(49) 



(50) 



The values of u m and v m depend on details on the Zee- 
man coupling and the electron-electron interaction and 
can be fixed by energy minimization. 

The excitations \S) , \AS) have collective coherence be- 
tween single-particle states with different spins and an- 
gular momenta values that differ by one. The RSR of the 
field S_ (r) = S x (r) — iS y (r) for the Skyrmion state is 
given by 



S_(r) = (5|*Ur)* t (r)|5> 



(51) 



2^ K,m+1 ( r ) V™ ( r ) KnVn 
m=0 



where ^J. (r) is the field operator that creates an electron 
at r with spin a. The corresponding density of the two 
spin components in the RSR are then 



n t (r) = (S|*}(r)* t (r)|S) 



(52) 



rn— 



IVm( r )l \ v mY 



and 



nUr) = (5|*t( r) * (r )|5) 



(53) 



I Vo (l - )! 2 + I Vm+l ( r )! 2 \ U r. 



m=0 

The total density is of course 

n (r) = n-f (r) + n± (r) , (54) 
while the z— component of the spin field is given by 



S z (r) = ~ [nt (r) - n x (r)] 



(55) 



(Note that for a filled level, \S Z (r)| = \/Att( 2 .) 

Figure U shows the density and spin patterns for a spin 
Skyrmion with Q = 1 added to a filled Landau level. We 
have chosen for this figure the simple expression 



10 + A 



m + 10 



10 + A' 



(56) 



where A is the Skyrmion size. This expression gives 
a density and pseudospin pattern for the Skyrmion 
which is qualitatively close to that given by the energy 
minimization^. Note that all spin and pseudospin densi- 
ties in this figure and the subsequent ones in this paper 
are in units of l/27r£ 2 . This corresponds to the density of 
a filled Landau level (see Eqs. (|46|47p ). 



B. Orbital Skyrmion 

We now consider the case where all electrons have spin 
up and we try to make an orbital Skyrmion by flipping 
some orbital pseudospins from n = to n = 1. For sim- 
plicity, we assume full valley and spin polarization so that 
we can drop the layer and spin indices. Full valley polar- 
ization is expected at all odd integer filling factors when 
the interlayer potential is strong. The ground state at 
v = 1 is given by 



\gs) = n 4,™ io> • 



(57) 



m=0 



An orbital anti-Skyrmion state can be written as 



oo 

i- 45 ) = n [ u ™ c im + Vmc i. 



m+2 



m— — 1 



|0>. (58) 
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FIG. 4: (Color online) A charge Q = 1 Skyrmion for the 
parameters of Eq. (|56|) with A = 4. (a) Spin profile in x — y 
plane, (b) Total density n(r). 



We see that the angular momentum difference Am = 
m\ — m — — 2 because the lowest value of m in n = 1 is 
m = — 1. 

For the Skyrmion excitation with the same vorticity, 
we have three choices corresponding to p = —1,0,1 in 
the expression 



oo 

n[- 



u mC\ ,m+2 v Tn c Gm 



4, P io> 



(59) 



The RSR of the pseudospin field S_ (r) is given by 

oo 

S - W = - H h hm+2 ( r ) h 0,m (r) (60) 



m=Q 



while the densities in n = and n = 1 are given by 
no(r) = (^|*S(r)* (r)|^,) 



(61) 



2 

m I i 



m=0 



and 



n 1)P (r) = ^^(r)*!^)^) (62) 

OO 

= |01,p Wl 2 + \ h l,™+2 Wl 2 l«m| 2 ■ 



m=0 



For the orbital Skyrmion, the iota? density is given by 



n (r) 



£ (,S p |*J(r)^(r)|5 p ) 

i,j=0 



(63) 



= n (r) + ni, p (r) + 2Rc [5_ (r) 
and includes the extra contribution 2Re [5_ (r) 



For the z— component 
1 



S z (r) = - [n (r) 



(64) 



To give an example, we take p = — 1 and choose again 
Eq. (|56l) with A = 0.05. The total density and pseu- 
dospin patterns are represented in Fig. [5] These pat- 
terns differ markedly from those of a spin Skyrmion. The 
coupling nil ~ mo — ~2 between the angular momenta 
in ii = and n = 1 makes the orbital-pseudospin vec- 
tor to rotate by 47T instead of 27r around the Skyrmion 
center. Because the value of u m is small, the profile of 
the ^-component of S z (r) is basically the (inverted) den- 
sity profile m (r) (see Fig. [3]). Also, the total density 
n (r) is anisotropic. Comparing this profile with that of 
N p (r) = no (r) + m tP (r) (not shown in the figure) which 
is isotropic, we understand that the anisotropy comes 
from the term 2Re [S- (r)] in Eq. (1551 . The density pro- 
file is typical of the density pattern of an electron in state 
n = 0,p = — 1 (see Fig. [3J. 

Our choice for the u m , v m in this example is motivated 
by the fact that it reproduces qualitatively the density 
we observed in the crystal phases that we discuss later. 
Of course, the correct values of u m , v m must be obtained 
by energy minimization i.e. by solving the eigenvalue 
equation for the Hartree-Fock Hamiltonian of an isolated 
Skyrmion. This is discussed in Ref. For the same 
reason, we choose in Eq. (|59p a pairing with Am = — 2 
because the vorticity of the Skyrmions we get in our crys- 
tals is 47T. 

It is clear from Eq. (|59l) that there are many varia- 
tions of the microscopic wavefunction \S P ) that we can 
make that would lead to pseudospin textures with dif- 
ferent topological and real electric charges. A study of 
these different solutions is, however, beyond the scope of 
this paper. 
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FIG. 5: (Color online) An orbital Skyrmion with p = — 1 
added to a filled Landau level for the parameters of Eq. (|56[) 
with A = 0.05. (a) Orbital pseudospin profile in the RSR. 
(b) Total density n (r) in the RSR. 



C. Skyrmion energy 

Orbital Skyrmions could occur, for example, at v = —3 
when the bias is strong enough for all the charge to be 
transferred to the state \K', 0) . When d = 0, this state 
occurs^ for As > 0.00T2e 2 /K£ In this state, the gap be- 
tween the n = and n — 1 state, /3As, is very small and 
one would expect orbital Skyrmions, in analogy with spin 
Skyrmions, to have lower energy than the correspond- 
ing electron quasiparticles. The energetics of an orbital 



Skyrmion is however different from its spin analog as we 
now show. 

With all electrons in the state \K',Q), the energy per 
electron is given by 



E 1 1 

N~o = -2 As -2 X °' OAo(0) 



(65) 



The energy needed to remove one electron in orbital m 
from the ground state (i.e. the energy to create one hole) 
is given by 



Ehole — + ^o,o : o,o (0) , 



(66) 



while the energy required to add one electron in n = 1 is 
given by 

E e = ~A B +PA B -Xo,i,i,o(0) (67) 

and is negative at zero bias. The energy to create an 
electron (in n = 1) and hole (with n = 0) pair with 
infinite separation is thus 



1 R 

eh, (orbital) = E e + Ehole = f3A B + ~ J — ■ 



(68) 



Let's compare this result with the energy needed to create 
an electron (in n = 0) and hole (in n = 0) pair with 
different spins i.e. 



eh, (spin 



9VbB 



(69) 



The first term on the right-hand side of this last equa- 
tion is the Zeeman energy. The energy needed to flip 
an orbital pseudospin is smaller than the energy to flip 
a spin (at zero bias and Zeeman couplings). It follows 
that the condition required to excite a Skyrmion pair i.e. 

A-skyrmion—antiskyrmion ^ ^eh, (orbital) IS more restric- 
tive for an orbital Skyrmion than for a spin Skyrmion. 
Periorbital) < &eh,(s P in) because of the presence of the 
extra exchange energy, Xo,i,i,o (0) i between different or- 
bitals which is not present in the spin case. 

Another important difference between the orbital and 
spin Skyrmions is that the exchange energy is smaller in 
n = 1 than in n = i.e. J ljUjl (0) = |^o,o,o,o (0) • 
For this reason, the gain in exchange energy obtained by 
making an orbital pseudospin texture is not as big as for 
a spin texture. 

A detailed numerical calculation of the energy of the 
Skyrmion and anti-Skyrmion excitations under finite bias 
thus has to be made in order to compare their energy 
with those of the electron and hole excitations. This can 
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by done by energy minimization, using the method de- 
scribed in Ref. 5> Despite our efforts, we have not been 
able so far to achieve sufficient precision with our numeri- 
cal code to classify energetically these different solutions. 
We thus concentrate, in this paper, on crystal solutions 
which are a lot easier to compute. 

We remark that orbital Skyrmions have been stud- 
ied previously in a conventional semiconductor 2DEG 16 . 
They were called, in this context, inter-Landau-level 
Skyrmions and involved spin flips between the n = 
spin down state and the n — 1 spin up state. Since 
there is no exchange energy between states with differ- 
ent spin indices, the energetics of these inter-Landau-level 
Skyrmions is different from that of our orbital Skyrmions 
in which orbital pseudospin flips occur between states 
with the same spin. The conclusion of Ref. [ID that 
inter-Landau-level Skyrmions are never the lowest lying 
charged excitation cannot be applied to our system. 

VI. SKYRMION CRYSTALS 

We now consider the ground state of the 2DEG in bi- 
layer graphene at non integer filling factors v € [1,3]. 
We look for crystal solutions of Eq. (|25|) allowing for the 
possibility of both valley and orbital pseudospin textures. 
We do not attempt to make an exhaustive study of the 
phase diagram of the 2DEG since considering all the pos- 
sible crystal states would be a formidable task. We re- 
strict ourselves to square and triangular lattices with one 
and two electrons per unit cell and single out the state 
with the lowest energy. We take £i = and discuss the 
effect of a finite £i in the conclusion of this paper. 

All results presented in this paper are for a magnetic 
field B = 10 T. We measure the energy in units of 
e 2 / 'k£ = 0.036 eV for k = 5 appropriate for graphene on 
a Si02 substrate. For the validity of the two-band model 
used in our calculations, we need Ab « 71 = 0.39 eV 
i.e. A s / (e 2 /n£) « 11. 

We start by presenting the crystal state at large bias 
for v around 1 because the Skyrmion at each site is close 
to the simple solution we presented in Fig. [S] in this case. 
We then study the crystals for v near 1 and 3. These two 
filling factors give very similar solutions. We end with 
filling factor near 2 where radically different solutions 
are obtained. 



A. Orbital Skyrmion crystal at large bias 

The simplest crystalline structure occurs at large bias 
with v around 1. In this case, valley K is empty and 
the charge is entirely in valley K' . This corresponds to 
the situation we studied in Sec. V(b). The crystal solu- 
tion is a triangular lattice of orbital Skyrmions with one 
Skyrmion per lattice site. We show an example of this so- 
lution for v = 1.2, Ab/ (e 2 /nl) = 1.28 in Fig. El We see 
that for each Skyrmion in the lattice, the pseudospin and 



density profiles (in the RSR) are close to the very crude 
Skyrmion solution we illustrated in Fig. [5j Note that if 
we use the GCR instead, the pseudospin and density pro- 
file for each Skyrmion are exactly those of the usual spin 
Skyrmion we illustrated in Fig. [4] i.e. the pseudospins 
rotates by 27r around the center of the Skyrmion and the 
pseudospins point downward at the center (we show this 
in Fig. E|) . 



(a) 




FIG. 6: (Color online) Orbital Skyrmion crystal at v = 1.2 
and Ab/ (e 2 /k£) = 1.28. (a) Pseudospin texture in the RSR. 
(b) Total density n (r) in the RSR. 

The Wigner crystal solution (i.e. no orbital pseudospin 
texture) can be found if the bias is taken to be extremely 
large i.e. of the order of Ab/ (e 2 / k£) ~ 30 which is 
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well beyond the limit of validity of our model. If we 
compare the interaction energy of the Skyrme crystal at 
As/ (e 2 /n£) — 1.28 with that of the Wigner crystal at 
Ab/ (e 2 //-c£) = 30, we find that the energy of the former 
is lower than that of the later by approximately 0.5%. 
(The interaction energy includes all terms in Eq. (|28|) 
with the exception of the bias energy). The Hartree 
part of the total energy is bigger in the Skyrme than 
in the Wigner crystal while the exchange (Fock) energy 
is more negative in the Skyrme crystal probably because 
Skyrmions are larger objects and overlap more with their 
neighbors. It could be also that isolated Skyrmions have 
lower energy than electron or hole quasiparticles. As 
mentioned earlier, we have not been able to confirm that 
Skyrmions have lower energy than isolated electrons and 
holes in the dilute limit using separate isolated quasipar- 
ticle calculations. If indeed Skrymions are only stable 
beyond a minimum density, the crystal stability must be 
related to inter- Skyrmion exchange energies 



The crystal state at v = 3.2 and zero bias is identi- 
cal to that at v = 1.2 with the exception that it is the 
electrons in the \AS, 1) that now condense into a crystal 
phase and orbital Skyrmions are formed by flipping or- 
bital pseudospins from \ AS, 0) to \ AS, 1). There are again 
2 orbital Skyrmions per site. The two states \S, 0) , \S, 1) 
are completely filled and inert and give a background 
density of 2/2ir£ 2 . The pseudospin and density patterns 
in the RSR and GSR for this crystal state with lattice 
spacing a are shown in Fig. [7J (In these figures, a is the 
Skyrmion lattice constant and a » do). The pseudospin 
pattern for Sk,z ( n ot shown in the figure) is identical to 
that of Sk'.z. The additional central peak in the RSR 
density profile occurs because there are 2 Skyrmions per 
site in this crystal. Note that the pseudospin for the 
charge q = 2e Skyrmion rotates by 4ir in the RSR but 
only 2?r in the GCR. 



Bias A B > A ( g } 



B. Skyrme crystals near v — 1,3 and zero bias 

It was shown in Ref. [l] that, in the Hartree-Fock ap- 
proximation, the ground states of the 2DEG in N = at 
integer filling factors satisfy a set of Hund's rules in which 
the spin polarization is maximized first, then the layer 
polarization is maximized to the greatest extent possi- 
ble, and finally the orbital polarization is maximized to 
the extent allowed by the first two rules. At zero bias, 
the ordering of the first four states (with spin up) is given 
by 



|5,0) 
15,1) 
\AS, 0) 
\AS,1) 



(70) 



71^ 1) + 71 

' \K,0) ' 



V2 
1 

71 



\K,1) 



V2 
1 

71 



K\ I) . 
K',0). 
K', 1) . 



in this order. (The guiding-center index X is left implicit 
in these equations). 

At v = 1, the first level is completely filled while at 
v = 3, the first three levels are completely filled. At 
v = 1 + x with |ar| < 0.5, a finite density of electrons 
(x > 0) or holes (x < 0), in the \S,l) or \S, 0) state 
condense into a crystal phase. At x = 0.2 and zero bias, 
we find that the electrons in the \S, 1) state condense into 
a triangular crystal with two orbital Skyrmions per site 
(i.e. one orbital Skyrmions in each layer). There is no 
layer-pseudospin texture (no pseudospin rotation) in this 
case since all electrons are in a symmetric state of the 
bilayer but there is an inter layer coherence. Moreover, 
an orbital-pseudospin texture is always present. 



With finite bias, the charge in layer K is progressively 
transferred to layer K' . For the crystal states discussed 
above, that means that the size of the Skyrmions de- 
creases in layer K and increases in layer K' . In our 
mean-field approximation, the charge of the Skyrmions 
is not quantized and the crystal states can be seen as ex- 
otic charge density waves with complex pseudospin tex- 
tures. Above a very small bias of order / (e 2 / ntj sa 
0.0011 at v = 1, all electrons are pushed into the 
\K',0) valley and interlayer coherence is lost. The or- 
dering of the energy levels at v = 1 is then given by 
\K', 0) , \K', 1) , \K, 0) , \K, 1) and the ground state has all 
electrons in \K', 0). At v = 1.2 and A B / (e 2 /n£) = 0.002, 
we find a triangular crystal of orbital Skyrmions in layer 
K' with again two electrons per site. There are no elec- 
trons in valley K. The guiding-center density and vector 
fields patterns for this state are given in Fig. [HI The total 
density is identical to that of the crystal state found at 
zero bias (see Fig. [7]) but it is now completely in layer K' 
instead of being equally shared between the two layers. 
The orbital pseudospin pattern in K' is the sum of the 
orbital pseudospin patterns found in each layer at zero 
bias. 

The pseudospin texture and density for a single orbital 
Skyrmions of charge q — 2e on each site of the lattice in 
Fig. [7] can be obtained with the microscopic expression 



5(2e) ) = II [-^ C U+2 + «m4, ro ] c\ fl c\,-l |0) 
m=0 

The density in this state is given by 



r4 2e) (r) = (r)| 2 + |4> 1)0 (r) 



(71) 
(72) 



\ h hrn+2 (r)| 2 \V. 



2 i |2 




FIG. 7: (Color online) Orbital Skyrmion crystal at v — 3.2 and As/ (e 2 /«£) = 0. (a) Orbital pseudospin texture in the RSR. 
(b) Orbital pseudospin texture in the GCR. (c) Total density n (r) in the RSR. The pseudospin pattern for Sk,z (not shown) 
is identical to that of Sk',z- (d) Total density n (r) in the GCR. The lattice constant is a. 



while the pseudospin texture is still given by Eq. (1601) and 



Si 2e) (r) 



n (r) 



,(2e) 



(r) 



(73) 



The phase diagram of the liquid state at v = 3 is 
richer than at v = 1 since a mixed state with both or- 
bital and interlayer coherences is possible 3 due to the fact 
that the kinetic energy contribution — PoAb in Eq. © 
is negative and so decreases the exchange-enhanced gap 
A* = E n= o < E n= i with increasing bias. (This gap A* 
is positive at zero bias but changes sign at sufficiently 
high bias). In the region with orbital coherence only, 



i.e. for As/ (e 2 / 'k£) > 0.0022, the ordering of the levels 
is given by \K' , 0) , \K', 1) , \K, B) , \K, AB) where B and 
AB represent bonding and anti-bonding combinations of 
the n = and n — 1 states defined by 



\K,B) 
\K,AB) 



with 



,(2)' 



(74) 



^\K,0)+VT^\K,1), (75) 



(76) 
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(a) 




FIG. 8: (Color online) Orbital Skyrmion crystal at v = 1.2 
and As/ (e 2 /k£) = 0.002. (a) Orbital pseudospin texture in 
the RSR. (b) Total density n (r) in the RSR. The pseudospin 
S K = 0. 



and 

(at = 10 T) is the bias at which the electrons are com- 
pletely transferred to the state \K, 1) at v = 3. The two 
levels \K',0) , \K' , 1) are completely filled. The number 
of electrons in state \K, 1) is given by 



VK A = -^y (78) 

At v = 3.2 and A B > A^ c) (with = 0.0021, we 

get a triangular crystal of orbital Skyrmions with charge 
q = — 2e per site with a density and pseudospin pat- 
terns in the n = 0, 1 basis similar to those represented 
in Fig. [U The only difference with the v = 1.2 case is 
that Skyrmions are now made from a filled \K, B) level 
by flipping pseudospins to the \K, AB) level. When A# 

is close to Ag', however, electrons are in majority in 
level n = and there is not much difference with the 
v = 1.2 case. When the bias is sufficiently strong for the 
exchange-enhanced gap to be negative, the crystal phase 
is much more complex but the orbital pseudospin tex- 
ture persists. We do not discuss this limit further in this 
paper. 

Skyrmion with charge q = — 2e are not unknown. It 
was shown in Ref. [13, for example, that at small den- 
sity, there is an attractive force between two Skyrmions 
with opposite global phases of their spin component that 
goes like 1/R where R is the separation between the 
two Skyrmions. At large separation R, this force pre- 
vails over Coulomb repulsion. Also, in previous studies 
of spin and pseudospin Skyrmions in conventional semi- 
conductor's 2DEG, it was found that lattice with pairs 
of Skyrmions occurred for small value of the Zeeman or 
bias couplings 12 . 



D. Skyrme crystals near v = 2 

At v = 2, the uniform ground state at zero bias has 
intcrlayer coherence in n — and in n = 1 so that the 
first two states in Eqs. (1701) are filled. Above a critical 

bias of the order of A^ c) / (e 2 / ' k£) w 0.003, all charge 
are transferred to the K' valley and states n = and 
ii = 1 are then fully occupied. Interlayer as well as orbital 
coherences are lost. 

At v = 2.2 and zero bias, the electronic phase consists 
of two layer-pseudospin meron crystals with each meron 
carrying charge q = — e/2. There is a meron texture in 
Po and in Pi. This phase is depicted in Fig. [9j The 
merons are arranged in a checkerboard configuration with 
8 merons per unit cell (the total charge in one unit cell 
is 4e). The layer-pseudospins in the central meron are 
rotated by a phase ir with respect to the merons at the 
corner of the unit cell. The orbital coherence is more 
than ten time smaller than the interlayer coherence and 
can be neglected. 

The charge in the pseudospin merons with pseudospin 
down (up) at the center progressively decreases (in- 
creases) when the bias is increased. For v — 2.2 and 
at As / (e 2 / K £) w 0.007, we find that there is a phase 
transition to a state where the two states in valley K' 
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FIG. 9: (Color online) Interlayer Skyrmion crystal at v — 2.2 
and As/ (e 2 /k£) = 0.(a) Layer-pseudospin texture in n — 
and (b) Layer-pseudospin texture n = 1 in the RSR. (c) Total 
density n (r) in the RSR. 



are completely filled and the remaining electrons crystal- 
lize in the K valley. In the n = 0, 1 basis, this gives a 
triangular crystal with one electron per site and an or- 
bital pseudospin vortex around each electron. Note that 
the liquid phase at v = 2 has all electrons in a bonding 



state \K,B) of Eq. ([74]) (with = {v - 1)^/tt/2/4P 
in this case) and we could have expected a Wigner crys- 
tal phase with electrons in the \K, B) state at each site. 
It seems however that the system again prefers to form 
a pseudospin texture at each site. This crystal state is 
represented in Fig. [TU] 



(a) 




x/a 

(b) 




FIG. 10: (Color online) Orbital crystal at v = 2.2 and 
As/ (e 2 /ni) = 0.007. (a) Orbital pseudospi n texture in the 
RSR. (b) Total density n (r) in the RSR. 

As in the v — 3.2 case discussed above, the exchange- 
enhanced gap A* between n — and n — 1 is of the 
order of e 2 / 'n£ at v — 2. The contribution — (3qAb (see 
Eq. |(8} ) decreases this gap as is increased. This 
increases the number of flipped orbital pseudospins. At 
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A^ = (y— l)-y/7r/2/4/3, there is a transition to a Wigner 
crystal with all electrons in n — 1 at each site and there 
is no pseudospin texture anymore. 

It is possible to find another interesting solution with 
our numerical code at Ab > that has, however, an 
higher energy than that of Fig. [lOj We mention it here 
because it is closely related to the work reported in Ref. @. 
This solution is a crystal of layer-pseudospin Skyrmions 
in n = and in n = 1 with a charge q — — 2e Skyrmion at 
each crystal site. This solution is the natural extension 
of the solution at zero bias since the bias transfers the 
charge of half the merons in one layer to the other half in 
the same layer. This structure is shown in Fig. QTJ Wc 
remark that our convention for the interlayer pseudospin 
is that state up corresponds to valley K. Since the charge 
is pushed in valley K' with positive bias, the majority 
state is pseudospin down. In that case, a Skyrmion has 
spin up at the center and spin down away from the center. 

This type of solution i.e. Skyrmions with a superpo- 
sition of n = and n = 1 interlayer textures have been 
studied by Abanin et al£ (the small contribution /3qAb 
to the gap was set to zero in that paper). These authors 
concluded that such charge q = —2e Skyrmions would 
have lower energy than electron or hole quasiparticles at 
filling factor v — 2.0. The crystal structure that we get 
is consistent with their finding but it is not the ground 
state. As we mentioned before, our conclusions for the 
crystal state do not necessarily apply to the case of an 
isolated Skyrmion. If we define a Skyrmion creation op- 
erator in states n = 0, 1 as 
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m=—l 



(80) 



then the 2e Skyrmion state can be written as 

SM)=44\0) (81) 



and the angular momentum pairing is such that mx 1 — 
itik = 1 for both Skyrmions. 

In a previous publication^, we derived an effective 
model for the orbital pseudospin-wave excitations at 
v = 3. This effective model had in it a Dzyaloshinskii- 
Moriya interaction. This type of interaction favors the 
formation of spiral or vortex states. We see that this 
term is also at work in the crystal states. 
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FIG. 11: (Color online) Interlayer-Skyrmion crystal at v = 
2.2 and Ab/ (e 2 /k£) = 0.03. (a) Layer-pseudospin texture in 
n = and (b) n = 1 in the RSR. (c) Total density n (r) in the 
RSR. Note that the majority state is pseudospin down. This 
state is not the ground state in our numerical calculation. 



VII. TOTAL AND LOCAL DENSITY OF 
STATES 



Skyrmion lattices with charge q — — 2e can be distin- 
guished from Skyrmion lattices with charge q — — e by 
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their total density of states (TDOS) which is defined by 
1 

7T 



9T M 



7T 



/ rfrIm[G^ a (r,r, W ) (82) 



where G 



(i?)a,a 



(with a = K,K' and n = 0, 1) is the re- 
tarded single-particle Green's function which is related 
to the Matsubara Green's function defined in Eq. (f2"3")l 

by (<1> cu + iS) = G i n i i a,a (q, w) . 



(a) 

40 



30 



(/> 
'E 

20 

CO 

o 
a 

H 10 



-2.0 



-1.8 

E/(e 2 /Kl) 



-1.6 




E/(e 2 /Kl) 



FIG. 12: Total density of states for orbital Skyrmion crystals 



at v = 1.2. (a) For A s /(e 2 /K 



1.28. Skyrmions with 



q = -e. (b) For A B / (e 2 / 'ni) = 0.002. Skyrmions with q = 
— 2e. Only the low-energy part of the TDOS is shown in these 
figures. 

The number of peaks near the Fermi level in the TDOS 
is equal to the number of electrons in a Skyrmion. This is 
illustrated in Fig. [T^] where we show the low-energy part 
of gr (w) corresponding to the crystals of Fig. [6] (q = — e 
Skyrmion at each site) and Fig. [JJ (q — — 2e Skyrmion 
at each site). A similar result was also found for bubble 
crystals in semiconductor's 2DEG^&. 



It was shown by Poplavskyy et alJ£ that the density 
pattern in the bubble crystal can also be seen by scanning 
tunneling microscopy. This measure is related to the local 
density of states (LDOS) which is defined by 



9l (r,w) 



-i^Im (r,r, W )" 



(83) 



ttS* 



where 



G a n % (q,«) 



dre lq - r G^>' a (r,r,w), (84) 
N^ n (-q,oj)K n , n (q). 



We show in Fig. [13] the LDOS in valley K 1 evaluated at 
the energy of the two highest-energy peaks in Fig. [T2T a) 
and at the highest-energy peak in Fig.[T2Tb). The LDOS 
is almost the same for both peaks in the case of the 
Skyrmion crystal with charge q — — 2e. The LDOS for 
the Skyrmion crystal with charge q = — e looks much 
the same. Following Ref. 113, we can also sum the LDOS 
evaluated at all the peaks below the Fermi energy. It 
is easy to show analytically, using Eq. ( |24|) . that this 
summation gives 



E F 



g L (r, w) du = N p (r) , 



(85) 



where N p (r) is the density we defined in Sec. V(b). This 
density is actually quite close to the real space density 
n (r) that we plotted in many of the figures of this paper 
(it does not contain the term 2Re [S- (r)]). 



VIII. ELECTRIC DIPOLE TEXTURES 

Apart from the minus sign in front of {p a ,y (G)) in 
Eq. (|4"Tj) , the vector field representation for the electric 
dipoles in the crystal states with orbital coherence is 
just like the GCR of the orbital pseudospin field S a ± (r) 
where a = K, K' . We give an example of the dipole field 
in Fig. [7] for the charge q = — 2e orbital Skyrmion crys- 
tal at v = 3.2 and bias Ab/ (e 2 /Kt) = 0. The rotation 
of the pseudospins is 2tt for both charge q = — e and 
charge q = —2e Skyrmions so that the S a ,j^ field pattern 
does not allow to discriminate between these two types 
of crystals. 

From Eq. (|4"0|) , we see that in the presence of an exter- 
nal uniform electric field E = Eo x x + Eo y y in the plane 
of the layers, the coupling with the electron gas is given 
by 

H ext = V2£eN v J dr (i:,, J> h , (r) - E 0tV p j>y (r)) (86) 
= -V2leN v (Eo, xPj , x (G = 0) — E 0>yPj>y (G = 0)) . 
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IX. CONCLUSION 
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FIG. 13: (Color online) Local density of state (LDOS) in val- 
ley K' for the orbital Skyrmion crystals considered in Fig. [12] 
(a) and (b). (a) and (b) LDOS for the two higher-energy 
peaks of the crystal with 2 electrons per site, (c) LDOS at 
the higher-energy peak for the crystal with one electron per 
site. 



In the liquid phase where the orientation of the orbital 
pseudospins in the x — y plane is arbitrary, this term 
allows us to rotate the orbital pseudospins in that plane. 
For a Skyrmion crystal, the effect may be more complex. 
The parallel electric field forces the orbital pseudospin in 
the x — y plane and so should increase orbital coherence. 
It should also change the form of the orbital pseudospin 
texture of the Skyrmions, orienting them more towards 
the x axis. 



We have presented in this paper a study of some crys- 
tal phases with valley and/or orbital pseudospin textures 
that can occur in bilayer graphene away from integer fill- 
ing factors in Landau level N = 0. Our calculations are, 
strictly speaking, valid within the two-band tight-binding 
model introduced in Sec. II and within the Hartree-Fock 
approximation. 

In our numerical calculatons, we have neglected the 
terms 04 and A in Eq. (Q]). These terms change the 
gap between the two orbital states n — and n = 1 
to Ci — PoA-b in the K valley and to (i + (3qAb in the 
K' valley (see Fig. [5]). As we mentioned in Sec. II, the 
value of £i is not known precisely. If we take the values for 
Po, (3 4 and A cited in Sec. II, we find Ci/ (e 2 / ' ni) = 0.113 
at B = 10 T so (i is probably not small. Since the critical 
bias needed to push the charge in one layer is such that 
PoA^ /{e 2 /K£) w 0.177 x 10 -3 , we see that these addi- 
tional terms have the possibility to change the phase dia- 
gram in an important way especially the phases at small 
or zero bias. Furthermore, the orbital coherence depends 
on the gap between the two orbital state. With the value 
of Ci cited above, the bias A# needed to place the orbital 
n=l below n = in the valley K is Ab « 1.27 i.e. a 
large value. Fortunately, our numerical calculations show 
that the orbital Skyrmions crystal (which is the most im- 
portant state we discussed in this paper) does survive in 
the phase diagram even with a finite £i ■ At filling factor 
v = 1.2, the additional gap suppress the q = — 2e orbital 
Skyrmion crystal in favor of a q — — e orbital Skyrmion 
crystal. This is consistent with our mention in Sec. VI 
(b) that Skyrmion with q = — 2e are found at small gap. 
At filling factor v = 3.2, the gap Ci — PoAb can be made 
small (or even negative) and orbital Skyrmion crystals of 
both types are found. 

Skyrmion crystals have both phonon and spin (or pseu- 
dospin) wave modes. In Ref. |2(X it was shown that the 
classical (or quantum mean-field) energy of the Skyrmion 
is independent of the angle ip which defines the global XY 
orientation of the spin components. This extra U (1) de- 
gree of freedom for a single Skyrmion leads to a broken 
symmetry in the crystal ground state and hence to a 
spin wave mode which remains gapless in the presence 
of a Zeeman field. We expect a similar gapless mode 
for a crystal of orbital Skyrmions. Fluctuations due to 
phonons and to these gapless pseudospin modes will have 
to be considered at finite temperature in order to evalu- 
ate the stability of the crystal structures discussed in the 
present paper. 
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